Advanced Regressions

Panel Data

There are a few packages out there that implement panel data models. We’ll cover two of them, plm() - an older frequently used package, and a newer package, fixest(). Another package to note is lfe.

Fixed Effects Models

Fixed effects models are equivalent to including dummy variables for the entity that you are following over time. You can include these dummy variables as factors (either for time or entity). However, that is computationally inefficient if you have some 1,500 people, for example. So, instead, we use the within transformation.

In plm() is an option for “within” - that is the fixed effects model. This is because you are taking differences within an entity.

Let’s use world bank data to apply the fixed effects model.

library("WDI")
library("tidyverse")

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──

## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.6     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.1.1     ✓ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

wdi <- WDI(country = "all", start=2000, end=2015, extra="TRUE",
           indicator=c("NY.GDP.MKTP.KD.ZG","SL.TLF.CACT.FE.ZS","SP.POP.TOTL", "SE.PRM.CUAT.ZS" ))

Let’s quickly clean this up a bit by renaming the columns. I’m going to leave in unnecessary columns since the data is not very large, but if I was working with a large dataframe, I’d drop it.

wdi <- rename(wdi,  gdp= NY.GDP.MKTP.KD.ZG, lfpart_f= SL.TLF.CACT.FE.ZS, pop= SP.POP.TOTL, edu = SE.PRM.CUAT.ZS)

names(wdi)

##  [1] "iso2c"     "country"   "year"      "gdp"       "lfpart_f"  "pop"      
##  [7] "edu"       "iso3c"     "region"    "capital"   "longitude" "latitude" 
## [13] "income"    "lending"

Let’s great a baseline model regression with OLS:

ols <- lm(gdp~edu + log(pop), data=wdi)

Now, let’s run our model.

We need to tell plm what entity is being followed and what is the variable for time using the option index, and we need to tell plm what kind of model we are running.

Typically, the function is a class of regression methods and there will be an option for a specific estimation model. plm can be considered a package for panel data models and the estimation model is “within”

A one-way fixed effects model is like this:

library("plm")

## 
## Attaching package: 'plm'

## The following objects are masked from 'package:dplyr':
## 
##     between, lag, lead

fixed <- plm(gdp~  edu + log(pop) , data=wdi, index = c("iso2c"), model="within")

summary(fixed)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = gdp ~ edu + log(pop), data = wdi, model = "within", 
##     index = c("iso2c"))
## 
## Unbalanced Panel: n = 143, T = 1-16, N = 625
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -17.802655  -1.053216   0.029157   1.430126  18.143716 
## 
## Coefficients:
##           Estimate Std. Error t-value Pr(>|t|)    
## edu      -0.044892   0.046240 -0.9708 0.332115    
## log(pop) -9.154496   2.422198 -3.7794 0.000177 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    6632.7
## Residual Sum of Squares: 6305
## R-Squared:      0.049406
## Adj. R-Squared: -0.23577
## F-statistic: 12.4738 on 2 and 480 DF, p-value: 5.2334e-06

#equivalent in stata:
#xtset iso2c year
#xtreg gdp edu lnpop, fe

But, we can include fixed effects for time, too. We just need to include it in the index.

fixed <- plm(gdp~ edu + log(pop) , data=wdi, index = c("iso2c", "year"), model="within")
summary(fixed)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = gdp ~ edu + log(pop), data = wdi, model = "within", 
##     index = c("iso2c", "year"))
## 
## Unbalanced Panel: n = 143, T = 1-16, N = 625
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -17.802655  -1.053216   0.029157   1.430126  18.143716 
## 
## Coefficients:
##           Estimate Std. Error t-value Pr(>|t|)    
## edu      -0.044892   0.046240 -0.9708 0.332115    
## log(pop) -9.154496   2.422198 -3.7794 0.000177 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    6632.7
## Residual Sum of Squares: 6305
## R-Squared:      0.049406
## Adj. R-Squared: -0.23577
## F-statistic: 12.4738 on 2 and 480 DF, p-value: 5.2334e-06

We can see the number of observations (N), the number of groups (n), how many units you have over time (T).

If you want to see the fixedeffects estimtates, you can do this with fixef()

fixef(fixed)

##      AD      AE      AL      AM      AO      AR      AU      AW      AZ      BA 
## 112.529 148.618 145.534 146.651 163.424 158.820 162.405 111.076 157.886 143.062 
##      BB      BD      BE      BF      BG      BH      BI      BN      BO      BR 
## 122.863 180.604 153.734 156.210 153.692 132.562 152.141 125.784 154.247 180.771 
##      BS      BT      BY      BZ      CD      CI      CL      CM      CO      CR 
## 123.381 129.265 151.712 121.774 176.574 165.433 161.109 159.106 169.218 148.298 
##      CU      CV      CY      CZ      DE      DK      DM      DO      DZ      EC 
## 154.809 123.908 132.740 154.571 172.647 146.775 106.038 156.089 163.745 158.763 
##      EE      EG      ES      ET      FJ      FR      GE      GH      GM      GN 
## 144.079 175.286 166.588 179.614 127.702 170.127 148.554 166.663 138.996 152.982 
##      GR      GT      GY      HK      HN      HR      HU      ID      IL      IN 
## 153.208 156.102 127.939 151.934 152.789 145.416 155.630 185.591 153.282 199.306 
##      IQ      IT      JM      JO      JP      KE      KG      KH      KR      KW 
## 169.719 168.001 142.170 152.737 179.463 168.857 149.168 158.552 171.036 143.675 
##      KY      KZ      LB      LI      LK      LS      LT      LU      LV      MD 
## 106.757 161.270 153.603 102.857 155.697 140.120 144.389 132.375 139.265 143.961 
##      ME      MH      MK      ML      MN      MO      MR      MT      MU      MV 
## 129.221 103.889 138.063 157.516 143.076 136.348 144.256 129.504 137.036 144.114 
##      MX      MY      MZ      NA      NE      NG      NL      NO      NP      OM 
## 175.439 165.288 163.939 135.399 164.001 180.609 157.648 148.010 161.910 142.607 
##      PA      PE      PH      PK      PL      PR      PS      PT      PW      PY 
## 145.306 166.537 177.464 179.713 167.141 144.147 149.070 152.254  90.532 150.089 
##      QA      RO      RS      RU      RW      SA      SC      SE      SG      SI 
## 147.401 161.172 150.979 181.043 157.105 164.166 109.209 153.708 150.740 138.654 
##      SK      SL      SN      SR      ST      SV      SY      TD      TG      TH 
## 149.089 149.555 152.659 130.020 118.287 147.738 161.312 167.515 151.508 172.832 
##      TJ      TO      TR      TT      TZ      UA      UG      US      UY      VE 
## 155.887 116.845 175.294 132.738 168.767 175.605 166.069 185.291 145.614 165.566 
##      WS      ZA      ZW 
## 120.784 168.695 156.563

fixest

You can run the same above regression with the package fixest using the function feols(). I’m mentioning this particular package because it is insanely fast. It can handle a huge amount of fixed effects where stata might break (unless you use reghdfe). It can support non-linear models, high-dimensional fixed effects, multiway clustering and a bunch of options that come in handy when you work with detailed data or big data. So, I want you to know it exists because of it’s flexibility in all things fixed effects.

Let’s try it. Rather than setting the fixed effects with an index, you’ll include the variables that are fixed after |

Let’s see what I mean here:

library("fixest")

fixed_est = feols(gdp~  edu + log(pop) | country + year, data = wdi)  ## Fixed effect(s) go after the "|"

## NOTE: 3,631 observations removed because of NA values (LHS: 249, RHS: 3,623).

fixed_est

## OLS estimation, Dep. Var.: gdp
## Observations: 625 
## Fixed-effects: country: 143,  year: 16
## Standard-errors: Clustered (country) 
##           Estimate Std. Error   t value Pr(>|t|) 
## edu      -0.010863   0.051449 -0.211140  0.83308 
## log(pop) -4.105289   3.085700 -1.330424  0.18551 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 2.65713     Adj. R2: 0.476812
##                 Within R2: 0.008502

Standard errors are already clustered by country automatically, but if you want standard errors, you can do so using summary and the option se

summary(fixed_est, se = 'standard')

## OLS estimation, Dep. Var.: gdp
## Observations: 625 
## Fixed-effects: country: 143,  year: 16
## Standard-errors: IID 
##           Estimate Std. Error   t value Pr(>|t|)    
## edu      -0.010863   0.043952 -0.247154 0.804898    
## log(pop) -4.105289   2.262458 -1.814526 0.070241 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 2.65713     Adj. R2: 0.476812
##                 Within R2: 0.008502

There might be time where you want to cluster your standard errors with specific items. You can also specify what you would like to cluster with the option cluster

summary(fixed_est, cluster = c('iso2c'))

## OLS estimation, Dep. Var.: gdp
## Observations: 625 
## Fixed-effects: country: 143,  year: 16
## Standard-errors: Clustered (iso2c) 
##           Estimate Std. Error   t value Pr(>|t|) 
## edu      -0.010863   0.051449 -0.211140  0.83308 
## log(pop) -4.105289   3.085700 -1.330424  0.18551 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 2.65713     Adj. R2: 0.476812
##                 Within R2: 0.008502

Random Effects Models

For plm, changing to a random effect model is pretty simple. You just need to change the option “model” to “random”

re <- plm(gdp~edu + log(pop), data=wdi, index = c("iso2c"), model="random")
summary(re)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = gdp ~ edu + log(pop), data = wdi, model = "random", 
##     index = c("iso2c"))
## 
## Unbalanced Panel: n = 143, T = 1-16, N = 625
## 
## Effects:
##                  var std.dev share
## idiosyncratic 13.136   3.624 0.799
## individual     3.305   1.818 0.201
## theta:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1061  0.2940  0.4465  0.3877  0.4848  0.5539 
## 
## Residuals:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -17.9553  -1.7391  -0.0018  -0.0137   1.6605  23.4716 
## 
## Coefficients:
##               Estimate Std. Error z-value  Pr(>|z|)    
## (Intercept) 11.6506226  2.3066628  5.0509 4.398e-07 ***
## edu         -0.0587729  0.0099989 -5.8779 4.155e-09 ***
## log(pop)    -0.1870426  0.1257400 -1.4875    0.1369    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    9343.1
## Residual Sum of Squares: 8436.8
## R-Squared:      0.097062
## Adj. R-Squared: 0.094159
## Chisq: 34.6357 on 2 DF, p-value: 3.0126e-08

There is something else called a mixed model where you might have both fixed effects and random effect parameters. This is particularly useful for nested datasets and when you might want to include fixed effects at one level and random effects for another variable. I’ve rarely seen it in economics papers, but it is often included in psychology models. I won’t discuss this indepthly, but here is a website that explains how to apply it in R - although there are many out there and often people use the package lmer().

Hausman Taylor test

You may know that when deciding between random or fixed effects models, you have to ensure that the random variables do not correlate with your independent variables. Frequently, this assumption is broken and is why most people use fixed effects models. However, we often want to test between our fixed or random effect model. We do this with the phtest() function in plm.

fe <- plm(gdp ~ edu +log(pop) , data=wdi, index = c("iso2c"), model="within")

summary(fe)

## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = gdp ~ edu + log(pop), data = wdi, model = "within", 
##     index = c("iso2c"))
## 
## Unbalanced Panel: n = 143, T = 1-16, N = 625
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -17.802655  -1.053216   0.029157   1.430126  18.143716 
## 
## Coefficients:
##           Estimate Std. Error t-value Pr(>|t|)    
## edu      -0.044892   0.046240 -0.9708 0.332115    
## log(pop) -9.154496   2.422198 -3.7794 0.000177 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    6632.7
## Residual Sum of Squares: 6305
## R-Squared:      0.049406
## Adj. R-Squared: -0.23577
## F-statistic: 12.4738 on 2 and 480 DF, p-value: 5.2334e-06

phtest(fe, re)

## 
##  Hausman Test
## 
## data:  gdp ~ edu + log(pop)
## chisq = 17.116, df = 2, p-value = 0.000192
## alternative hypothesis: one model is inconsistent

There are a variety of tests for serial correlation, heteroskedasticity, random effects (if you should include random effects compared to an ols model, stationarity/unit roots, etc). These are outlined in the vignettes of plm and you can explore them further on your own.

Logit/Probit

Logit/probit models can be run using the generalized linear model packages glm() You identify a logit or probit model using the option “family”. You need to specify the link function to distinguish between a logit or probit model.

A good introduction to logit models (and many other applied econometric models) can be found at UCLA’s website here and for probit models here

Let’s start with a logit model:

setwd("/Users/mkaltenberg/Documents/GitHub/Data_Analysis_Python_R/Advanced Regressions Stargazer/")
smoking <- read.csv("smoking.csv")
names(smoking)

##  [1] "smoker"   "smkban"   "age"      "hsdrop"   "hsgrad"   "colsome" 
##  [7] "colgrad"  "black"    "hispanic" "female"

logit <- glm(smoker ~ age + female + hsdrop +hsgrad +colsome +colgrad, data = smoking, family = "binomial")
summary(logit)

## 
## Call:
## glm(formula = smoker ~ age + female + hsdrop + hsgrad + colsome + 
##     colgrad, family = "binomial", data = smoking)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.0622  -0.8203  -0.5884  -0.4056   2.2899  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.940401   0.135308 -14.341  < 2e-16 ***
## age         -0.005998   0.001964  -3.054 0.002256 ** 
## female      -0.210057   0.048467  -4.334 1.46e-05 ***
## hsdrop       1.771164   0.126802  13.968  < 2e-16 ***
## hsgrad       1.521847   0.113672  13.388  < 2e-16 ***
## colsome      1.161500   0.116063  10.007  < 2e-16 ***
## colgrad      0.424145   0.125638   3.376 0.000736 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 11074  on 9999  degrees of freedom
## Residual deviance: 10583  on 9993  degrees of freedom
## AIC: 10597
## 
## Number of Fisher Scoring iterations: 4

The probit model just needs a tweak in the link function within the option “family”, and that is set like this:

probit <- glm(smoker ~ age + female + hsdrop +hsgrad +colsome +colgrad, data = smoking, family = binomial(link = "probit"))
summary(probit)

## 
## Call:
## glm(formula = smoker ~ age + female + hsdrop + hsgrad + colsome + 
##     colgrad, family = binomial(link = "probit"), data = smoking)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.0563  -0.8214  -0.5919  -0.4005   2.3070  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.125266   0.073137 -15.386  < 2e-16 ***
## age         -0.003458   0.001157  -2.990 0.002791 ** 
## female      -0.123336   0.028445  -4.336 1.45e-05 ***
## hsdrop       1.005005   0.068663  14.637  < 2e-16 ***
## hsgrad       0.851358   0.059047  14.418  < 2e-16 ***
## colsome      0.637198   0.060379  10.553  < 2e-16 ***
## colgrad      0.220119   0.064848   3.394 0.000688 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 11074  on 9999  degrees of freedom
## Residual deviance: 10584  on 9993  degrees of freedom
## AIC: 10598
## 
## Number of Fisher Scoring iterations: 4

Marginal effects can’t be read directly from the output. You need to apply the link function to interpret the results (though the sign and standard error can be interpreted as usual).

This is why the margins function is so important! Margins works with other regressions, as well.

Margins

The margins package is quite similar to the margins package in STATA. You can apply interpretation to a wide variety of regressions - logit/probit, but also non-linear terms in OLS models. You can check out more examples and features in its vignette here.

The package makes it easy to calculate values and plot them.

When we use margins(model_name) it will give us the average marginal effect for each variable

library(margins)

logit_m <-margins(logit)

And we can plot it quit easily:

plot(logit_m)

To get the partial effect at the mean (or at any particular value), we us the at option:

Here we specify a range of numbers

margins(logit, at = list(age = c(18,40,60)), variables = "age")

Or at the mean of age

mean_age = mean(smoking$age, na.rm=TRUE)

margins(logit, at = list(age = c(mean(smoking$age, na.rm=TRUE))), variables = "age")

Marignal effects for non-linear terms in OLS

The margins package can also be used to get the marginal effects for non-linear terms in OLS models - often these are interactions or polynomial functions. Let’s take an example of a model with the squared term population.

Notice that to square a term I use “^2” and the “I()” to indicate that it’s a second order term.

int <- lm(gdp ~ lfpart_f + pop + I(pop^2), data = wdi)
summary(int)

## 
## Call:
## lm(formula = gdp ~ lfpart_f + pop + I(pop^2), data = wdi)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -65.787  -2.051   0.003   2.061 119.426 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.281e+00  2.993e-01  10.964  < 2e-16 ***
## lfpart_f     1.250e-02  5.594e-03   2.234 0.025541 *  
## pop          1.289e-09  2.510e-10   5.134 2.98e-07 ***
## I(pop^2)    -1.836e-19  4.822e-20  -3.807 0.000143 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.109 on 3636 degrees of freedom
##   (616 observations deleted due to missingness)
## Multiple R-squared:  0.01023,    Adjusted R-squared:  0.009412 
## F-statistic: 12.53 on 3 and 3636 DF,  p-value: 3.807e-08

Let’s apply marginal effect to population. We can take the marginal effect of any particular value of population that we are interested in.

# in python wdi['pop']

m_pop <- mean(wdi$pop, na.rm=TRUE)  # getting the mean of population

margins(int, at = list(pop =m_pop))

We can combine those two lines of code into one (and not save the value of the mean of population that will take up RAM)

margins(int, at = list(pop = mean(wdi$pop, na.rm=TRUE)))

We can input a list of values rather than just one value

margins(int, at = list(pop = c(1000000,5000000, 9000000)))

It’s far better to use numbers relevant to the value/data that exists. We can use a bunch of summary statistics instead of just one to see the marginal effects of population. Turkey’s 5 numbers include min, lower hinge, median, upper hinge and maximum - the function is called fivenum()

Looking at the marginal effects, you can see that the labor force participation of women is constant (as expected), but population marginal effects depend on the value of population - there seems to be a non-linear relationship.

fivenum(wdi$pop, na.rm= TRUE)  # This is how to get the statistics

## [1]       9392    1336126    9088884   58056732 7339076654

margins(int, at = list(pop =fivenum(wdi$pop, na.rm= TRUE)))  #integrated within margins command

## Warning in check_values(data, at): A 'at' value for 'pop' is outside observed
## data range (97962,7339076654)!

And we can also graph the marginal effect of population with cplot()

cplot(int, "pop", what = "effect", main = "Average Marginal Effect of Population")

Stargazer (Again)

Let’s show some models side by side

library(stargazer)

## 
## Please cite as:

##  Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.

##  R package version 5.2.2. https://CRAN.R-project.org/package=stargazer

knitr::opts_chunk$set(echo = TRUE, warning = FALSE)

stargazer(ols, re, fe, type = "html", column.labels = c("OLS", "RE", "FE"), model.names = FALSE,
          dep.var.caption = "",
          title            = "Panel Data Results",
          covariate.labels = c("Edu", "lnPop", "Constant"),
          dep.var.labels   = "GDP per capita",
          out = "panel_results.html")

## 
## <table style="text-align:center"><caption><strong>Panel Data Results</strong></caption>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td colspan="3">GDP per capita</td></tr>
## <tr><td style="text-align:left"></td><td>OLS</td><td>RE</td><td>FE</td></tr>
## <tr><td style="text-align:left"></td><td>(1)</td><td>(2)</td><td>(3)</td></tr>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Edu</td><td>-0.059<sup>***</sup></td><td>-0.059<sup>***</sup></td><td>-0.045</td></tr>
## <tr><td style="text-align:left"></td><td>(0.008)</td><td>(0.010)</td><td>(0.046)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">lnPop</td><td>-0.275<sup>***</sup></td><td>-0.187</td><td>-9.154<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.093)</td><td>(0.126)</td><td>(2.422)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">Constant</td><td>13.052<sup>***</sup></td><td>11.651<sup>***</sup></td><td></td></tr>
## <tr><td style="text-align:left"></td><td>(1.747)</td><td>(2.307)</td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td></tr>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>625</td><td>625</td><td>625</td></tr>
## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.084</td><td>0.097</td><td>0.049</td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.081</td><td>0.094</td><td>-0.236</td></tr>
## <tr><td style="text-align:left">Residual Std. Error</td><td>4.082 (df = 622)</td><td></td><td></td></tr>
## <tr><td style="text-align:left">F Statistic</td><td>28.677<sup>***</sup> (df = 2; 622)</td><td>34.636<sup>***</sup></td><td>12.474<sup>***</sup> (df = 2; 480)</td></tr>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td colspan="3" style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>

Lecture 12: Advanced Regressions & Stargazer Continued

We’re going to cover: